############
##Libraries#
############
suppressMessages(library(plot.matrix)) #plot correlation matrix
suppressMessages(library(RColorBrewer)) #color palette
suppressMessages(library(R.matlab)) #import mat files
suppressMessages(library(dplyr)) #mean and sd by group
suppressMessages(library(splitstackshape)) #random stratified sampling
suppressMessages(library(reshape2)) #plotting libraries
suppressMessages(library(ggplot2))
suppressMessages(library(gridExtra))
suppressMessages(library(grid))
suppressMessages(library(egg))
suppressMessages(library(pbapply)) #progress bar in apply functions#Choose rat and session
rat <- "rat08"; folder <- "190410"; wind <- 128
########################
#Load all tetrode files#
########################
setwd(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/", rat, "/", rat, "_", folder,"_statistics", sep=""))
listFrames <- list.files(pattern="spikesFrame.*?\\.csv")
spikesTetr <- lapply(listFrames, read.csv)
#Get number of units detected in every tetrode
tetrUnits <- sapply(spikesTetr, function(x) length(unique(x$unit)))
#Change the unit number in the full data frame
sumIni <- c(0, cumsum(tetrUnits)[-length(tetrUnits)])
for (i in 1:length(tetrUnits)) {
spikesTetr[[i]]$unit <- spikesTetr[[i]]$unit + sumIni[i]
}
##########################
#Load noise removal times#
##########################
setwd("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal")
noiseRem <- readMat(paste(rat, "_", folder, "_window_", wind, ".mat", sep=""))[[1]]
#1=noise, 0=no noise
noiseThres <- floor(0.8*length(spikesTetr))
int2rem <- colSums(noiseRem) >= noiseThres#Function to label each spike as noise ot not
noiseMark <- function(events, noiseMat, wind) {
#1=noise, 0=no noise
int2rem <- colSums(noiseMat) >= noiseThres
#Time interval vector
timesInt <- (0:(dim(noiseMat)[2]+1))*wind
#Find to which interval corresponds each spike event
whichInt <- sapply(events$spikes*1e3, function(x) findInterval(x, timesInt))
#Label each spike as noise or not
noise <- sapply(whichInt, function(x) as.integer(int2rem[x]))
#Store in data frame
spikes2denoise <- data.frame(spikes = events$spikes, unit = events$unit, stimuli = events$stimuli, interval = whichInt, noise = noise)
spikes2denoise <- spikes2denoise[complete.cases(spikes2denoise), ]
return(spikes2denoise)
}
# Store all tetrodes in a list
list2denoise <- pblapply(spikesTetr, function(x) noiseMark(x, noiseRem, wind))
saveRDS(list2denoise,
file=paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal/", rat, "_", folder,
"_NoiseRemoval_window", wind, ".RData", sep=""))
# list2denoise <- readRDS(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/noise_removal/", rat, "_", folder,
# "_NoiseRemoval_window", wind, ".RData", sep=""))
#Check percentages
for (i in 1:dim(noiseRem)[1]) {
print(table(list2denoise[[i]]$noise))
}
0 1
2677 18017
0 1
7153 9687
0 1
1527 18275
0 1
14948 17581
0 1
1266 16883
0 1
7078 19308
0 1
7652 8106
0 1
8471 17496
tetrs <- 1:dim(noiseRem)[1] #Rat01, 02, 08
threshold <- c("1", "005", "01", "2", "1", "01", "05", "2") #Rat08 190410
# tetrs <- c(1:5, 8) #Rat 03
# threshold <- c("005", "1", "05", "01", "005", "4") #Rat03 190306
# threshold <- c("2", "01", "005", "015", "03", "005", "05") #Rat02 190306
# threshold <- c("05", "2", "025", "025", "10", "2", "8", "3") #Rat01 190308You must enable Javascript to view this page properly.
#Join all tetrodes into one data frame
spikesFull <- dplyr::bind_rows(list2denoise); spikesFull <- spikesFull[order(spikesFull$unit), ]
#Stimuli trains
trains <- readMat(paste("D:/Sergio/MSI Prestige/PhD/GCM_rats/IPHYS_David/data/", rat, "/Trains_", rat, "_", folder,
"_RAW.mat", sep=""))
feeder <- unlist(trains[[1]][1][[1]][2])
lever <- unlist(trains[[1]][1][[1]][3]); leverSec <- rep(c("press","release"), length(lever)/2)
#CAUTION: IF DATA IS "MERGED", SOFTWARE IS THE 5TH ELEMENT, IF NOT, 4TH
software <- unlist(trains[[1]][1][[1]][4])
#0 begin/end baseline, 1 reward, 2 non reward left, 3 non reward right, 4 end behavioral
stimOrder <- c(0,0,1,2,3,2,1,3,2,1,3,1,3,2,3,1,2,3,1,2,1,2,3,2,1,3,2,1,3,
1,3,2,3,1,2,3,1,2,1,2,3,2,1,3,2,1,3,1,3,2,3,1,2,3,1,2,4)
stimCount <- c(1:57)countsFinal <- list()
#Parameters
binSize <- 2
freq <- 1e4
spikesNoNoise <- spikesFull[spikesFull$noise == 0, ]
maxTime <- floor(max(spikesNoNoise$spikes)*freq)
for (g in 1:3) {
#Subset stimuli
spikeSubset <- spikesNoNoise[spikesNoNoise$stimuli == g, ]
minSpike <- min(spikeSubset$spikes); maxSpike <- max(spikeSubset$spikes)
#Generate full time vector with freq precision
vecTimes <- seq(0, maxTime, by = 1)
stim2 <- floor(software[stimOrder == g]*freq); stim2end <- stim2 + 30*freq
#Storage objects
countsToCorr <- list() #for spike counting
####COUNT SPIKES####
for (i in unique(spikesNoNoise$unit)) {
#Spikes time to freq precision
freqSpikes <- floor(spikeSubset[spikeSubset$unit == i, 1]*freq)
#Spike detection in full time vector
spikesPresence <- vecTimes %in% freqSpikes
##
stim2Count <- vector()
for (j in 1:length(stim2)) {
stim2Count <- c(stim2Count, spikesPresence[stim2[j]:stim2end[j]])
}
##
##
sumInt <- seq(1, length(stim2Count), binSize*freq)
binCounts <- vector()
for (k in sumInt[-length(sumInt)]) {
binCounts <- c(binCounts, sum(stim2Count[k:(k+(binSize*freq))]))
}
##
countsToCorr[[i]] <- binCounts
}
#Store in one matrix
lenMat <- length(countsToCorr[[1]])
countsToCorr.matrix <- matrix(unlist(countsToCorr), ncol = lenMat, byrow=T)
countsToCorr.matrix[is.na(countsToCorr.matrix)] <- 0
countsFinal[[g]] <- countsToCorr.matrix
}
uu <- countsFinal[[1]]; vv <- countsFinal[[2]]; yy <- countsFinal[[3]]; zz <- cbind(uu,vv, yy)
corMat <- cor(zz)
spear.corMat <- cor(zz, method=c("spearman"))
kend.corMat <- cor(zz, method=c("kendall"))####SPIKE COUNTS PLOT######
rastermat <- cbind(uu, vv, yy)
rastermat[rastermat > 30] <- 30 #adjust for every rat and session
melted_rastMat <- melt(rastermat)
midpoint <- (min(rastermat)+max(rastermat))/2
meanpoint <- mean(rastermat)
medianpoint <- median(rastermat)
ggplot(data = melted_rastMat, aes(x=Var2, y=Var1, fill=value)) +
ggtitle("Spike count") +
geom_tile() +
scale_fill_gradient2(low = "red", high = "yellow", mid = "white", midpoint = midpoint,
limit = c(min(rastermat),max(rastermat)), space = "Lab", name="Spikes count") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = (9))) +
scale_y_continuous(trans = "reverse") +
# geom_hline(yintercept=c(5.5,8.5,10.5,12.5,15.5,17.5), linetype="solid", color = "green", size=1) +
geom_vline(xintercept=c(dim(uu)[2]+0.5, 2*dim(uu)[2]+0.5), linetype="solid", color = "blue", size=2)colfunc.shape <- colorRampPalette(c("navy", "white", "red"))
numcols <- 100
colcol.shape <- colfunc.shape(numcols)
plot(rep(1,numcols),col=(colcol.shape), pch=19,cex=2)nunits <- dim(uu)[1]
#Format stim 1 mat
binSize <- 2
numtot <- 18*(30/binSize)
sepVec <- seq(1, numtot + 5, 30/binSize); sepEnd <- sepVec[-1] - 1
sepList <- list()
for (i in 1:length(sepEnd)) {
sepList[[i]] <- uu[, sepVec[i]:sepEnd[i]]
}
uu.RealStim <- do.call(rbind, sepList)
uu.Neurons <- rep(1:nunits, 18)
uu.toRM <- data.frame(unit = factor(uu.Neurons), uu.RealStim)
# aa <- rmcorr(unit, "X1", "X2", uu.toRM)
#Format stim 2 mat
sepList <- list()
for (i in 1:length(sepEnd)) {
sepList[[i]] <- vv[, sepVec[i]:sepEnd[i]]
}
vv.RealStim <- do.call(rbind, sepList)
vv.Neurons <- rep(1:nunits, 18)
vv.toRM <- data.frame(unit = factor(vv.Neurons), vv.RealStim)
#Format stim 3 mat
sepList <- list()
for (i in 1:length(sepEnd)) {
sepList[[i]] <- yy[, sepVec[i]:sepEnd[i]]
}
yy.RealStim <- do.call(rbind, sepList)
yy.Neurons <- rep(1:nunits, 18)
yy.toRM <- data.frame(unit = factor(yy.Neurons), yy.RealStim)
#################uu.toRMmeans <- matrix(unlist(by(uu.toRM[,-1], uu.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)
vv.toRMmeans <- matrix(unlist(by(vv.toRM[,-1], vv.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)
yy.toRMmeans <- matrix(unlist(by(yy.toRM[,-1], yy.toRM$unit, colMeans)), nunits, (30/binSize), byrow = TRUE)
zz.toRMmeans <- cbind(uu.toRMmeans, vv.toRMmeans, yy.toRMmeans)
corMat.toRMmeans <- cor(zz.toRMmeans, method="kendall")
meancolDiagKen <- corMat.toRMmeans
diag(meancolDiagKen) <- 0library(reshape2)
library(ggplot2)
melted_cormat <- melt(meancolDiagKen)
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) +
geom_tile() +
scale_fill_gradient2(low = "navy", high = "red", mid = "white", midpoint = (min(meancolDiagKen)+max(meancolDiagKen))/2,
limit = c(min(meancolDiagKen),max(meancolDiagKen)), space = "Lab", name="Kendall\nCorrelation") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1)) +
coord_fixed() +
geom_hline(yintercept=c((30/binSize)+0.5, (60/binSize)+0.5), linetype="solid", color = "green", size=1) +
geom_vline(xintercept=c((30/binSize)+0.5, (60/binSize)+0.5), linetype="solid", color = "green", size=1)#Comparison labels
division1 <- dim(uu.toRMmeans)[2]
division2 <- dim(uu.toRMmeans)[2] + dim(vv.toRMmeans)[2]
division3 <- dim(uu.toRMmeans)[2] + dim(vv.toRMmeans)[2] + dim(yy.toRMmeans)[2];
labs <- c(rep(c(rep("s1s1", dim(uu.toRMmeans)[2]), rep("s1s2", dim(uu.toRMmeans)[2]), rep("s1s3", dim(yy.toRMmeans)[2])), division1),
rep(c(rep("s2s1", dim(vv.toRMmeans)[2]), rep("s2s2", dim(vv.toRMmeans)[2]), rep("s2s3", dim(yy.toRMmeans)[2])), division1),
rep(c(rep("s3s1", dim(vv.toRMmeans)[2]), rep("s2s3", dim(vv.toRMmeans)[2]), rep("s3s3", dim(yy.toRMmeans)[2])), division1))
labMat <- matrix(labs, division3, division3, byrow = FALSE)
halfLabs <- labMat[lower.tri(labMat, diag=FALSE)]
halfCorMat <- corMat.toRMmeans[lower.tri(corMat.toRMmeans, diag=FALSE)]
#Mean, SE and SE
fullComp.means <- by(halfCorMat, halfLabs, mean)
fullComp.sd <- by(halfCorMat, halfLabs, sd)
fullComp.se <- by(halfCorMat, halfLabs, function(x) sd(x)/sqrt(length(x)))
#CI95
alpha <- 0.05
t <- qt((1-alpha)/2 + .5, 108-1)
fullComp.CI95 <- t*fullComp.se
#All in data frame
fullComp.frame <- data.frame(Comp = unique(halfLabs), Mean = as.vector(fullComp.means),
SD = as.vector(fullComp.sd), SE = as.vector(fullComp.se), CI95 = as.vector(fullComp.CI95))
#Barplot
library(ggplot2)
ggplot(fullComp.frame) +
geom_bar( aes(x=Comp, y=Mean), stat="identity", fill="forestgreen", alpha=0.5, width = 0.5) +
geom_errorbar( aes(x=Comp, ymin=Mean-SE, ymax=Mean+SE), width=0.2, colour="orange", alpha=0.9, size=1) +
ggtitle("using SE") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"))# dataStat <- data.frame(Comp = labs,
# Corrs = as.vector(corMat[, 1:division1]))
library(DescTools)
dataStat <- data.frame(Comp = halfLabs, Corrs = FisherZ(halfCorMat))
## parametric
summary(aov(Corrs ~ Comp, dataStat)) Df Sum Sq Mean Sq F value Pr(>F)
Comp 5 85.80 17.160 362.2 <2e-16 ***
Residuals 984 46.62 0.047
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.t.test(dataStat$Corrs, dataStat$Comp, p.adj = "bonf")
Pairwise comparisons using t tests with pooled SD
data: dataStat$Corrs and dataStat$Comp
s1s1 s1s2 s1s3 s2s2 s2s3
s1s2 < 2e-16 - - - -
s1s3 < 2e-16 1.00 - - -
s2s2 1.00 < 2e-16 < 2e-16 - -
s2s3 0.88 < 2e-16 < 2e-16 0.61 -
s3s3 6.2e-07 < 2e-16 < 2e-16 2.9e-07 8.5e-05
P value adjustment method: bonferroni
## non parametric
kruskal.test(Corrs ~ Comp, dataStat)
Kruskal-Wallis rank sum test
data: Corrs by Comp
Kruskal-Wallis chi-squared = 643.12, df = 5, p-value < 2.2e-16
pairwise.wilcox.test(dataStat$Corrs, dataStat$Comp, p.adj = "bonf")
Pairwise comparisons using Wilcoxon rank sum test
data: dataStat$Corrs and dataStat$Comp
s1s1 s1s2 s1s3 s2s2 s2s3
s1s2 < 2e-16 - - - -
s1s3 < 2e-16 1.00000 - - -
s2s2 1.00000 < 2e-16 < 2e-16 - -
s2s3 1.00000 < 2e-16 < 2e-16 1.00000 -
s3s3 0.00034 < 2e-16 < 2e-16 0.00043 0.00177
P value adjustment method: bonferroni
## stratified random sampling
stratDat <- stratified(dataStat, "Comp", 1000)
stratDat %>% group_by(Comp) %>% summarise(comp_mean = mean(Corrs), comp_sd = sd(Corrs))# A tibble: 6 x 3
Comp comp_mean comp_sd
<fct> <dbl> <dbl>
1 s1s1 0.956 0.257
2 s1s2 0.436 0.179
3 s1s3 0.412 0.168
4 s2s2 0.952 0.269
5 s2s3 1.01 0.249
6 s3s3 1.12 0.215
kruskal.test(Corrs ~ Comp, stratDat)
Kruskal-Wallis rank sum test
data: Corrs by Comp
Kruskal-Wallis chi-squared = 643.12, df = 5, p-value < 2.2e-16
pairwise.wilcox.test(stratDat$Corrs, stratDat$Comp, p.adj = "bonf")
Pairwise comparisons using Wilcoxon rank sum test
data: stratDat$Corrs and stratDat$Comp
s1s1 s1s2 s1s3 s2s2 s2s3
s1s2 < 2e-16 - - - -
s1s3 < 2e-16 1.00000 - - -
s2s2 1.00000 < 2e-16 < 2e-16 - -
s2s3 1.00000 < 2e-16 < 2e-16 1.00000 -
s3s3 0.00034 < 2e-16 < 2e-16 0.00043 0.00177
P value adjustment method: bonferroni